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Abstract 

The pair separation model of Goto and Vassilicos (S Goto and J C Vassilicos, 2004, New J.Phys., 
6, p. 65) is revisited and placed on a sound mathematical foundation. A DNS of two dimensional 
homogeneous isotropic turbulence with an inverse energy cascade and a k^'""^^ power law is used 
to investigate properties of pair separation in two dimensional turbulence. A special focus lies on 
the time asymmetry observed between forward and backward separation. Application of the present 
model to this data suffers from finite inertial range effects and thus, conditional averaging on scales 
rather than on time has been employed to obtain values for the Richardson constants and their ratio. 
The Richardson constants for the forward and backward case are found to be (1.066 ± 0.020) and 
(0.999 ± 0.007) respectively. The ratio of Richardson constants for the backwards and forwards case 
is therefore gb/gf ~ (0.92 ± 0.03), and hence exhibits a qualitatively different behavior from pair 
separation in three dimensional turbulence, where gt > p/ (J Berg et al. , 2006, Phys.Rev.E, 74(1), 
p. 016304). This indicates that previously proposed explanations for this time asymmetry based 
on the strain tensor eigenvalues are not sufficient to describe this phenomenon in two dimensional 
turbulence. We suggest an alternative qualitative explanation based on the time asymmetry related to 
the inverse versus forward energy cascade. In two dimensional turbulence, this asymmetry manifests 
itself in merging eddies due to the inverse cascade, leading to the observed ratio of Richardson 
constants. 

1 Introduction 

Lagrangian pair separation statistics are the vital ingredient necessary to calculate concentration fluc- 
tuations and concentration covariances which is important for studying pollution dispersal, combustion 
processes and the spread of contaminants in a liquid or gas [9, 13]. 

One of the fundamental quantities in turbulent particle dispersion is the mean square separation of 
an ensemble of particle pairs (A^(t)) [13], which will be studied in this article. 

Conversely, turbulent mixing, which is the convergence of particles that were further apart at an 
earlier time due to a turbulent flow, can be regarded as being equivalent to time inverted dispersion. 
Therefore, the underlying quantity of mixing processes is the so called backwards dispersion, that is the 
distribution of particle separations (A^(t)) at a time t for a prescribed separation Aq at a later time 
to > t [6, 10, 11, 21, 23]. 

Let us denote the forward mean square separation at timc^ t with an initial separation Aq at < = by 
(A^(t)]Ao)j^^. Equivalently, the backward mean square separation shall be denoted by (A^(— t)]Ao)|^^j. 
Note that the use of —t means that the time argument is actually positive, since for backward processes 
t < 0. For simplicity of notation in the backwards case, the minus sign shall be dropped henceforth and 
it is understood that time in the backwards case always refers to times smaller than the initial flow time 
t = 0. For flows whose dynamics are time- reversible, like Gaussian flows or kinematic simulations (KS), 
one could expect that the mean square separation of the forward and backward case coincide and indeed, 

^The time refers to the lapsed time of the underlying flow of the separation process, t = denotes the time when the 
initial / finial separation Aq is fixed. 



this has been shown for Gaussian flows and KS [10, 11, 21]: 

(A^WI^o),„,==(A^W|Ao),^,. (1) 

For other flows, like e.g. turbulent flows governed by the Navier-Stokes equation, this equality cannot 
be assumed and indeed, it has been shown using Lagrangian stochastic models [17, 21], experiment and 
DNS [3] that in three dimensional turbulence, backwards dispersion happens at a faster rate than forward 
dispersion. 

In this article we extend and give a sound mathematical foundation for the pair separation model 
as introduced by Goto & Vassilicos 2004 [14. shorthand GV04]. Wc start with a short summary of the 
GV04 model and then go on to refine and extend this model and explore its features and similarities to 
Richardson's distance neighbor function [20] . Finally, wc investigate the asymmetry between "forward" 
and "backward" pair separation in a DNS of two dimensional homogeneous isotropic turbulence, and 
compare the findings to those in three dimensional turbulence. 

2 Pair separation model 

The notion of pair separation as a process of burst-like separation events has been discussed and related 
to the streamline topology by Goto and Vassilicos [14] using a 2D DNS with an inverse 5/3 energy 
cascade. Wc shall remind the reader of the basic concept and introduce some refinements. 

In two dimensional multiscale flows, the basic idea is that each particle in a fluid is on a so called 
"patron" eddy. These patron eddies are coherent structures which persist for a time that is long enough 
to influence the dynamics of the fluid. They exist at every scale present in the flow and can be visualized 
by applying a coarse graining filter on the velocity field with an associated cut-off scale rjc- Each such 
eddy has typically an elliptic zero-acceleration point at its center. It is so called because the velocity 
gradient tensor du/ dx takes a vortical form around such a zero-acceleration point when one moves in the 
frame where the zero-acceleration point is at rest. Between the patron eddies, one finds hyperbolic zero- 
acceleration points which are locally surrounded by a straining velocity gradient tensor. Both elliptic 
and hyperbolic zcro-accclcration points arc assigned an associated scale rjc corresponding to the filtering 
scale. In the following, this scale rjc will be referred to as the "size" of the vortex or zero-acceleration 
point. 

The GV04 model portrays the pair separation process as a series of "bursts" to larger scales: Two 
distinct particles in a fluid share at least one patron eddy at each time (e.g. they always belong to the 
patron with the associated scale of rjc = where Lq is the size of the system / box size / etc.). Now 
consider the common patron eddy with the smallest associated scale r^c = A of a pair of particles at a 
particular time t. As long as both particles belong to this patron, their separation will also be typically of 
the scale A. It is now possible that at some time t H-Tj(A), one of the particles encounters a hyperbolic 
zero-acceleration point of scale A which removes it from the smallest common patron, where Ti^{-) is 
some function depending on the parameter ^. After such an encounter, the smallest shared patron of the 
particle pair will be of size ^A where ^ > 1. Hence, the typical pair separation will also have increased. 

While this picture strictly holds when locally moving with the zero-acceleration points and in two 
dimensions, it is possible to extend this notion to a more general context. Given sufficient persistence 
of the streamlines, GV04 argue that the presented arguments hold in a global reference frame when 
considering the elliptic and hyperbolic velocity stagnation points (i.e. points where the fluid velocity 
M = 0) instead of the zero-acceleration points in the local framc(s). This picture is not as intuitive as 
the patron notion in the local frames, however it seems like the natural generalization to a global frame. 
Hence, when considering sufficiently persistent multi-scale flows, also in three dimensions, one should be 
able to assume that particles belong to structures (previously called patrons) of a certain scale rjc- 

In three dimensions, these structures will not be flat eddies (like the patrons in two dimensions) 
but rather some kind of elongated eddies such as vorticity tubes or a patch of straining region in the 
velocity field. Despite not knowing the exact shape and properties of such persistent structures, they 
have frequently been observed in turbulence experiments and DNS [15, 22]. Therefore, one can assume 
that they are present and have a typical scale rjc which does not change much during the lifetime of such 
a persistent structure. 
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2.1 Refined model 



In the present work, we shall add to the model the possibility of two particles converging to a smaller 
characteristic separation A. While, as it turns out, there seems to be no practical benefit from this 
addition, it strengthens the derivation and mathematical foundation of the concept. 

2.1.1 Notation and basic ideas 

Let A„ = ^"Ao be the separation of a particle pair with initial separation Aq after a succession of 
separating and converging burst events, where n is an integer. Such a burst event can occur when 
at least one of the pair particles encounters a straining (hyperbolic) stagnation point. Note that ^ is 
the characteristic ratio of successive separations following a burst event, and as such is a constant that 
represents the respective ratio for burst events of all individual pairs. 

The probability 6„ of encountering a straining stagnation point of size A„ must be inversely pro- 
portional to the characteristic time T^(A„) between burst events for pairs with a separation w A„, as 
introduced by GV04 [14], and therefore proportional to the mean distance between hyperbolic stagnation 
points e{An) = ns{An)-^^'^: 

6„ocT^(A„)-i ^ 6„ (xu'n«(A„)i/^ (2) 
where the stagnation point density ns(A„) is given by [8] 

n,{fj,) = CsC-^ (^^) ' , (3) 

and C is the largest scale in a multi-scale flow (typically the integral scale), r]c is the coarse-graining 
scale and cannot be taken smaller than the smallest scale i] of the flow below which the power-law energy 
scaling fails (for this section's purposes, r]c = A„), d the dimension of the flow, Ds the fractal dimension 
of the stagnation point distribution in space and Cs is the stagnation point number which is proportional 
to the number of >C-sized stagnation points per unit volume. We assume a power spectrum E{k) oc k^f , 
1 < yO < 2 for (3) to be valid [8], in which case Ds has been shown to obey the relation [8] 

= (4) 

For the purposes of this paper, we assume p = 5/3. Due to the dependence on the stagnation point 
structure, one finds that the probability of encounter, 6„, mainly depends on the exponent of the energy 
spectrum p: 

for p = 5/3, and 

B = CbCI''^u'/C^'^ . (6) 

Cb is a proportionality constant for the probability of a particle pair to encounter a stagnation point. 

When a particle pair with separation A„ encounters a straining stagnation point of size A„, it is 
assumed to separate further with a probability p < 1, resulting in a separation of scale A„_|_i (separating 
burst). Alternatively, the pair can remain at the same scale of separation with probability 1 — p. This 
is the burst process described in GV04 [14], although the probability p was effectively absorbed into the 
coefficient B of (5) and not discussed by them. 

The converging burst process can be initiated when a particle pair with separation A„ encounters 
a straining stagnation point which is of smaller scale than the separation of the pair, i.e. A„_i. As 
mentioned previously in §2, both straining and elliptical stagnation points have an associated coarse- 
graining scale rjc- Between the eddies or structures of scale rjc = A„, and their elliptical stagnation points 
lie hyperbolic straining stagnation points of the same scale, which can act as a "gateway" into or out 
of a coherent structure of equal scale as shown in fig. 1. Hence, a particle pair of typical separation A„ 
needs to encounter a straining stagnation point of scale A„_i to experience the converging burst process 
bringing them together to a separation of A„_i. 

Hence, a particle pair with separation A„ encounters a straining stagnation point of scale A„_i with 
probability 5„_i, and then can converge to a separation of scale A„_i with a probability called q. 
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coherent structure of scale A„_i 

Figure 1: Diagram illustrating the pair convergence process. Shown is a particle pair and their tra- 
jectories. Their separation before encountering the stagnation point is A„. Reversing all arrows gives 
the figure for the separation process. Note that the scale of the stagnation point and of the associated 
coherent structure is always the same, no matter whether one considers the converging or separating 
process. 



2.1.2 Derivation of PDF of pair separation 

Analogously to GV04 [14]. the PDF of pair separation can be derived as follows. 

Let Qnit) be the probability for a particle pair to have a separation between A„ and A„-|_i at a 
certain time t and n is any integer. Qn will then evolve according to the probabilities 6„, 6„_i, p and q 
as given by the following evolution equation, 

^Qn = pbn-l Qn-l ^ {pK + <?6n-l) Qn+qbn Qn+1 , (7) 

at 

which is illustrated by fig. 2. Here, we have made use of the locality-in-scale hypothesis [13, 14] which 
suggests that an increase or decrease of Qn can only originate from the neighboring separation intervals 
associated with the probabilties Qn±i, but not from separation intervals that are further away, like e.g. 
the ones associated with Qn±i with i > 2. 




Figure 2: Diagram explaining the burstwise pair separation and convergence process. Note that if 
the arrows were reversed ("time-reversal"), the only resulting difference would be that p and q are 
interchanged. If an encounter with a hyperbolic stagnation point does not result in a burst or convergence 
event, the separation of a pair remains unchanged. 

In the limit of continuous separation A, we find that 

A„-.A(n)=rAo. (8) 

Defining a = In^, we arrive at 

Q„-^P(A„) = aAP(A), (9) 
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where f'(A) is the probability density function of continuous separation A. 

If one assumes that a is a small parameter it is possible to Taylor expand all Q,, and 6„ in (7), which 
results up to second order in a in 



1 dP 
B 'dt 



-a{p- q) + - a q 



_d_ 
OA 



Ai/3p 



p + q d 



2 dA\ dA 



A 



d 



Ai/3p 



(10) 



The index n has been dropped for legibility. This result reverts to the equation found by GV04 [14] when 
(J = and p = 1. This partial differential equation has 4 parameters, B, p, q and a. However, either p 
or q can be absorbed into B, leaving only 3 free parameters. In fact, the equation takes a tidier form if 
the rescaled time r is introduced, 

T ^^Ba^ {p + q)t, (11) 



along with the renormalized probabilities 



P = 



Therefore the equation now reads 



P 



p + q 



and 



p + q 



(12) 
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9 p* 
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d_ 
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/3 



P 



4 c^A I 9A ' 



(13) 



under the constraint p* +q* ~ 1, leaving effectively only two free parameters (e.g. q* and a). The third 
free parameter has been absorbed into the time variable, and is therefore implicitly still present. 

This form highlights the influence of the model parameters on the shape of the PDF: The two free 
parameters occur as a combination in the PDF, hence it will not be possible to determine both parameters 
from the shape of the PDF. This becomes even more obvious when the equation is written in the following 
form: 

dP , , .w. ap ^^^4/3g^ (14) 



^ = fcoA-2/3p- 

OT 



9A2 



where 



ki 
k2 



i p*-q 
2 a 

3 (fco + 1) ., 
9 

4 ' 



1 

4' 



(15) 
(16) 
(17) 



It is now clear that the solution P(A, r) in our model is solely determined by fcg = kQ{q* , a) for a given 
initial condition. The solution P{A,t), with time t rather than the scaled time r, is then given by an 
additional rescaling of the time axis. It is therefore apparent that even if it is possible to completely 
determine P(A, t) by experiment or DNS, this can only lead to the determination of two parameters. It 
will not be possible to determine all three model parameters a, q* and B from fitting data to this PDF. 



2.1.3 PDF of pair separation for arbitrary a 

While the assumption a 1 from the previous section leads to a compact formulation for the time 
evolution of the PDF of pair separation, it is a very restrictive assumption which basically states that 
the "burst" processes described earlier happen in a more or less continuous fashion, thus suggesting that 
there may be no bursts. 

We shall now present an argument that the PDF evolution of the form (14) is valid for all positive 
a, albeit the coefficients A:o,i,2 will be more complicated in terms of the model parameters a, p and q. 
First note that the sum over all probabilities is constant in time, 



(18) 



as can easily be verified from (7). Substituting the series with infinite number of terms 



Qn±l 



J=0 



A 



_d_ 
dA 



(AP) 



as well as (5) and (9) into (7) leads to the following evolution equation for the continuous PDF, 



B dt 

which can be cast into the form 



fe=i 



1 dP _ d 
B~dt ~ M 



OO 



.fe=0 



where the first coefficient is given by 



OO l. 

— a 



k=l 



kl 



p 



9 1 



and the remaining coefficients can be determined from the recursive relation 

k—n 



1 



OO i. 



k—n 



kl 



3=0 



2\3 

3 



k + 1 
n+1 



k + 1 
n + 1 



with 



^■n — 1 Y ^ 



being the Stirling numbers of the second kind. The first two coefficients, 

Co = 



9 3 4 

a a a 

a(p-q) + —{p + 5q) - —{p - 19q) + j-r^{p + Q5q) 
6 54 648 



a2 a3 ^4 

Ci = + (5p-llg) + — (21p + 85g) 



(19) 
(20) 



(AP) , (21) 



(22) 



(23) 



(24) 
(25) 

(26) 
(27) 



agree with (15) - (17) up to second order in a. Note that truncation after any k term in (22) preserves 
the property 

|/%(A,t) = 0, (28) 
which is the continuous-A form of (18). Furthermore, (22) can also be rearranged as follows 



1 dP 



B dt ~ ^ " 



m— 1 



A""2/^P(A) 



where the Dm are given by the recursive relation 

OO 

Dm = Cm-1 — E] Dk+1 
k—m 



k 

771—1 



with 



(a;)(") = x{x + l){x + 2) ■ • ■ (x + n - 1) 



{k-m+l) 



r{x + n) _ (x + n- 1)! 



r(^) 



(.T-l)! 



(29) 



(30) 



(31) 



6 



being the Pochhammer symbol or "rising factorial" . Again, the first two coefficients, 
D2 



-a{p - g) - y(3p - - ^(Sp -q)- ^(27p - 5g) 



-jip + + -ji'ip -q) + y^(2ip + 5<z) 



(32) 
(33) 



agree with (f5) - (17) up to second order in a and truncation after any m term in (29) preserves the 
property (28). 

One can now examine the time evolution of the second moment of P(A). Using (29) and two 
integrations by parts, one finds 



. 2 dP{^) 
dt 



dA 



poo poo 

-2 / i:>i A''/^PdA + 2 / DjA^'/^PdA 



9A"'-3 



Vu e {0,1,..., 00} 



(34) 

(35) 
(36) 



provided that all derivatives of P(A) involved in the third integral of (35) vanish sufficiently fast for 
A ^ 00: 

hm A-V3 ^ 

aA" 



A- 



0. 



(37) 



The conclusion to be drawn from this equation is that the complete evolution of (A^), and thus (A^) 
itself, depends only on the first two terms in (29) and therefore, can be determined from a Richardson- 
type equation such as (13). Hence, truncating the infinite series of (29) after the second term leads to 
an evolution equation for the PDF, which contains the zeroth, first and second derivatives of P{A) and 
gives the correct first and second moments only. Higher order moments of separation will require higher 
derivative terms to be included. Of course, the coefficients fco.1,2 of (14) will be functions of Di and D2, 
and thus of the model parameters B, p, q and a. Determining these 4 model parameters from the 3 
accessible coefficients fco,i,2 is generally not possible. 

Note that the coefficients fco.1,2 contain infinite series of m-terms a™/m!, and thus one can expect 
that these series will converge for arbitrary a. To avoid unnecessary complication of the notation, we 
shall use the fco,i.2 from the previous section, obtained under the assumption a ^ 1 for the remainder 
of this article. Note that this does not change the validity of the results, as they only depend on the 
finiteness of fco,i.2, and not on their dependence on the model parameters. 



2.1.4 Time reversal in the new model 

The caption of fig. 2 might hint that an exchange of the parameters p and q might be sufficient to account 
for a possible time asymmetry of pair separation. This, however, is misleading. 

The present model does not explicitly contain any dynamic features: all assumptions made can be 
satisfied by a Gaussian velocity field or a kinematic simulation, which is known not to exhibit the sought 
time asymmetry [10, 11, 21]. Hence, there is no reason to expect this model to account for the observed 
asymmetry. 

Instead, if one examines a velocity field which does exhibit the time asymmetry of mean square 
separation, it would imply that the model describes this situation as two separated cases with two sets of 
model parameters, one for forward separation and one for backward separation. It therefore can describe 
the asymmetry, but not explain it. 

Furthermore, one can consider the following mathematical argument. From experimental observations 
of the PDF of separation, one expects the solution P(A, t) to be "melting" , i.e. initial sharp peaks should 
flatten and spread out as time advances (in both the forward and backward case). This behavior is only 
observed when fco < 0. However (15) can only be negative if p* > q*. Hence, assuming that p* > q* in 
the forward case, fcp > in the backward case (or vice versa). 
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2.1.5 Comparison between the PDFs derived by Richardson, GV04 and the present work 

All three models [14, 20] can be represented by an equation of the form (14), where the time variable has 
been suitably normalized to satisfy fc2 = 9/4. Then, ki shows the same dependence on fco in all three 
models and the fcg are given by: 

present : fco = — 7: \- q + t (38) 

2 a 4 

GV04: ko = -^^ + 7 (39) 
2 aQY 4 

3 

Richardson: ko = -^(d-l) (40) 

The difference between the Richardson model and the other two is that for a given dimension of the 
problem (i.e. d = 2 or d = 3), the Richardson model predicts fep. Both other models do not predict the 
value of fco, but instead allow deviations from Richardson's PDF. GV04 have shown that their 2D DNS 
data is best fitted by acv ~ 1-3, whereas Richardson's PDF would predict acv = 1-5- 

Furthermore, for arbitrary a, we have shown in §2.1.3 that higher derivate corrections are necessary 
for higher separation moments. This could potentially explain deviations from Richardson's law for 
higher moments, should they be sufficiently resolved in experiments in the future. 

2.1.6 Richardson limit as guideUne for expected parameter values for q and a 

So far, data has been reasonably fitted with Richardson's form of the pair separation PDF [3, 14, 19]. 
One can therefore use the feg value in the Richardson case as a guideline for the expected parameters a 
and q (or q*; we shall use the ratio q/p instead of q* /p* purely for legibility). By equating (38) and (40), 
one can find a constraint on the parameters q and a to match the pair separation PDF as introduced by 
Richardson [14, 20]: 

9 _ _ l-a(l-i) 



p p* l + a(| + i) 



(41) 



One can invert (41) to find an expression for a: 



p{3d - 2) + q{M + 2) 
which leads to the following form for two and three dimensions 

3 p-q 



6__^Z:l___, (42) 



2 p + 2q 
6 p — q 
7p+^q 



for d = 2 (43) 
for d = 3 . (44) 



Note that for g = this matches the values found by GV04 [14]. 

2.1.7 Analytical solution for vanishing initial separation Aq = 

Equation (34) of GV04 [14], is an exact solution for (14), which we formulated here in terms of ko- It is 
valid for all three cases discussed here: 



-2ko 



P(A, r) - , (^] cxp 

^ ' ' r(3/2-2fco) V ^ 



(45) 



Note that this is only a valid solution when ko < as is it satisfied by e.g. Richardson and GV04. r(-) 
is the Gamma function. 

This analytical solution (see fig. 3) has the initial condition P(A, 0) = 25{A) (with 6{-) being the 
Dirac delta function) and hence, is not directly applicable to comparison with experiments as the (initial) 
separation between two physical particles is usually always > 0. However, if one assumes that after a 
certain time, the separation process will have "forgotten" the value of its initial state Aq, for r ^ (45) 
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(a) probability density P, ko = -3/4 





Figure 3: PDF of pair separation (45) with fco = -3/4 (a) and fco = -3/2 (b). Note that the PDF with 
ko = —3/2 spreads faster with respect to the scaled time r. See fig. 4 to compare the shape of both 
PDFs in A-direction. 
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can be compared to experiments with Aq > 0. This is frequently done in the hterature [3, 14, 19], and 
deviations from a fit are usually not discussed with sufficient care. 

This solution leads to the familiar t'^-law for the mean square separation, irrespective of the parameter 
fco, which only contributes to the coefficient Gr'- 



Jq 



(46) 



where 



(47) 



and 7 ~ — 2fco for brevity. More generally, one can find a formula for arbitrary moments of A with 
n> -(l + f7): 



(A"(r)) 



r 7- 



3 + 3 
2^2' 



r(7 + i) 



(48) 



Furthermore, normalizing the PDF of separation with its r.m.s. separation leads to a time-independent 
self-similar PDF [14]: 



where 



See fig. 4. 



P(A) = 



1/3 



3r(7 + 3/2) 



exp 



_^l/3^2/3 



P(A) = V(A2(t)} P(A, t) and A = A/V(A2(t)} 



(49) 
(50) 



probability density P 




Figure 4: Self-similar PDF P(A) for vanishing initial separation, fcg = —3/4 (solid) and ko = —3/2 
(dashed) is what Richardson [20] predicts respectively for 2D and 3D turbulence. 

The n-moments of this self-similar PDF are constants that only depend on 7: 

/ (A2)^ r(7 + f)^ r(7 + |)'"^ 

Provided one accepts the assumption a <C 1, one can expect these time- independent moments to be a 
good indicator of the validity of the presented model: the more moments we can fit to the data for a 
given 7, the closer this model represents what is going on. For arbitrary a, this relation is only valid for 
the first and second moment, as has been argued in §2.1.3. 

Since experimental data or data from numerical simulations cannot be obtained for Ao = 0, we need 
a means of understanding what happens for finite initial separations. 
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Ar2 






LUq 


dt 


V nil OL 7712 


30722 


680 


1.00147 


0.525 


6.1 X 10^5 


1.43 X lO--*"* 8 2.5 1 



Table 1: Parameters used in the numerical scheme as described in [14]: the number of grid points in the 
de-aliased scheme iV^, forcing wavenumber fc/, ratio of forcing wavenumber range /3, magnitude of fixed 
Fourier component of vorticity cJq, time increment of temporal integration dt, hyper- viscosity coefficient 
V and its exponent ttii, hyper-drag coefficient a and its exponent 7?i2 [14]. 

3 Two dimensional DNS study 

While the effect of time asymmetry of pair separation has been reported for three dimensional experiments 
and DNS [3, 17], there are no published results about this asymmetry in two dimensional turbulence. 

In this section, we shall present findings from the DNS run with a resolution oi N — 3072. We 
used the DNS scheme described in [14] with the parameters given in table 1 and typical scales given in 
table 2. Please refer to [14] for more details about this DNS, as reiterating the details would lengthen this 
article unnecessarily. All data were collected using 50000 particle pairs which were initially placed evenly 
distributed at random in the entire DNS domain. The distance between pair member particles is fixed 
at the distance Aq to give a delta function initial state of the pair separation PDF P(A, 0) = S{A — Aq). 
The spatial orientation of the separation vector was chosen at random. 



integral length C 




0.24 


smallest lengthscale in —5/3 intertial 


range rjf — 2-k /kf 


9.2 X 10-3 


outer vs. inner length scale L/rjf 




26 


r.m.s. velocity fluctuation u' ~ ^ {u ■ 


u)/2 


1.1 


integral timescale Tc — C/u' 




0.22 



Table 2: All quantities are given in DNS units and will be used to non-dimcnsionalizc where appropriate. 
Note that Tc ~ 3560 timesteps and hence, the DNS covers approximately 11 integral timescales. 

There are 20 runs in total: 10 for the initial separations, 

Ao e{^,|-,^,|-, 77/, 27?/, 477^,87;/, 1677^, 327?^} , (52) 

for both forward and backward separation in time. The particle trajectories are obtained by integrating 
the DNS velocity field u(f(t), t), 

x{t) ^xq+ I u{x{t'), t') dt' , (53) 
Jo 

using a second order predictor-corrector scheme [5] for 39150 consecutive time steps^ of the DNS which 
are stored on the hard drive. In the backwards case, the sequence of velocity field time frames is inverted 
and the velocity field negated to account for the reversed particle motion. 

3.1 Validity of model assumptions 

Before comparing the DNS data to the model presented earlier, it is necessary to check which part of the 
data actually satisfies the assumptions made in the model, namely the presence of a multiscale stagnation 
point topology which encompasses scales smaller and larger than the separation scale of the particle pairs 
under observation. 

Fig. 5 shows that the PDF of pair separation is rapidly spreading to larger separations, as would be 
expected from the exponential tail of the analytical solution (45). To quantify the time when many pairs 

^The number 39150 arose from hard disk constraints. 
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separation — 




Figure 5: Pair separation PDF P{A,t), plotted logarithmically in time. All figures show PDFs from 
forward separation data, the figures for backwards data look very similar, (a) to (c) have the following 
initial separations: Aq = rjf, 3277/. Black represents a PDF value of 10~^ or lower. The probability 
of a particle pair having a separation of the integral scale £ exceeds values of 10~^ (gray) at about 
t/Tc ~ 10-\ see text. 
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leave the inertial range, we introduce the threshold time tc which is defined as the time when P{C, t) 
crosses an arbitrary threshold for the first time: 



P{C,tc) = 10- 



(54) 



Note that the order of magnitude of tc and its dependency on Aq does not change significantly if a 
different threshold is chosen between 10~^ and 10~^. The threshold times for all given Aq show the 
expected decreasing trend in fig. 6. For comparison, the Batchelor time, 



1/3 



(55) 



indicating the end of the ballistic regime [3], is shown in the same figure. The dissipation rate e has been 
estimated from the constant energy flux of the DNS velocity field, as described in [14]. 
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Figure 6: Threshold time tc for all initial separations Aq. Forward separation data is shown as triangles 
(a), backwards data as diamonds (♦). For comparison, the Batchelor time is is represented by the 
dashed line. 

Summarizing, tc marks the time at which approximately 0.1%- 1% of pairs have a separation larger 
than the integral scale C and therefore have left the inertial range. The Batchelor time ts marks the 
time when initial separation effects are assumed to have ceased [2, 19] and after which the Richardson 
A oc t"^ law is expected. Also, separations below the forcing length rjf are outside the inertial range. 

Therefore, the presented model is strictly only valid for data with 

• rjf < Aq < £, i.e. initial separation within the inertial range, and 

• t < tc, i.e. enough pairs remain within the inertial range. 

Furthermore, the Richardson t'^ law can only be expected ior t > ts- In conclusion, the resolution of 
A'' = 3072 does not provide an inertial range that is wide enough to be able to observe the Richardson 
t^ law [3, 4, 24]. When comparing the data with the present model, it is therefore necessary to discuss 
the influence of these adverse effects in sufficient detail. 

The same applies when one wishes to compare the analytical solution (45) with no initial separation to 
DNS data: At the times the flt is performed, many pairs have left the inertial range already (^cutoff ^ tc) 
and hence, a basic model assumption is violated. To reflect this effect, the model would need to be 
extended by e.g. modelling the particle behavior similar to Brownian motion for pairs with A > C. 
However, this extension goes beyond the scope of this work, whose primary objective is to investigate 
the multiscale stagnation point topology within the inertial range. 



3.2 Mean square separation and parameter fit 

Despite the conclusion in §3.1 that it is not expected to observe a t^ law in the given DNS data, fig. 7 
shows a slope that is approximately close to t^ for times t G [ts, tc] for the runs with initial separation 
smaller than the forcing length: Aq < ??/. 
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Figure 7: Mean square separation from the DNS for all 10 initial separations. Forward separation is 
shown in black, backward in red. The Batchelor time ts is indicated by triangles (a) and the threshold 
time tc by diamonds (♦). The straight line is (x t^. The slope of the data is approximately oc in the 
intervals [tB,tc], although there is a noteable deviation for times just larger than ts- The 5**^ graph 
from the bottom has initial separation Ag = '7/- 

According to §3.1, this behavior was only to be expected for Aq > rjf. However, as mentioned previ- 
ously, the DNS resolution is currently not large enough to accommodate ts < tc for initial separations 
larger than rjf. 

Nonetheless, the approximate slope oc is quite meaningless, as for times t < tc, the mean square 
separation is still influenced by the finite initial separation. Thus, let us compare the model's prediction 
for (A^) with finite initial separation Aq > with the data and try to fit the model parameters. 

The model has been integrated for all 10 initial separations and for 3 different values of /cqj namely 
/sq e {—0.65, —0.75, —0.85}. The parameter r/t can be easily adjusted after integration of the model. 
This parameter has been manually adjusted to give a qualitative fit to the DNS data for the three values 
of fco mentioned previously. The two best fitting DNS data sets are shown in fig. 8 with Aq £ {r]f / A,rjf /2} 
for the best fit value of the time scahng parameter, r/t « 1.2. 

For initial separations in our data set which are smaller or larger than the best fit ones pointed out, 
the fit with these parameter values become less acceptable. Using different parameters improves the 
situation, but does not lead to a fit of the same quality as shown in fig. 8. 

Hence, we conclude that one set of parameters is not sufficient to fit all data for varying initial 
separations. 

From fig. 8, it is immediately apparent that the variation of fco only has a small impact on the variation 
of the curve and therefore, it is not possible to estimate fco to a very high accuracy from this kind of fit. 
Hence, while Richardson's value of fco is compatible with the present data, it can unfortunately not be 
uniquely confirmed. Nevertheless, by applying Ockham's razor, we shall assume Richardson's value of 
fco = —0.75 to be a reasonable fit for the remainder of this work. 

3.3 Dependence on initial separation 

The previous section highlights that the present model can not be fitted to the DNS data for all given 
initial separations. The validity of the model is thus clearly dependent on the initial separation, which 
is a further manifestation of the limited range of scales that the model operates on: In the DNS, too 
many particle pairs leave the range of multi-scale stagnation point distances too quickly and therefore, 
(A^) falls below the value expected from the present model. This effect is stronger for larger initial 
separations (see fig. 9). Hence, additional care needs to be applied when fitting data to the model, as 
the DNS data will deviate from the model for large times, as seen in fig. 9. 

For comparison, a simple functional form for the mean square separation with finite initial separation 
has been previously suggested to be used as a fitting law [14, 19]: 

(A^) = {[G^eY/H + AT)' = {g\'^t + ^.T)' . (56) 
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Figure 8: Best fit of numerically integrated mean square separation (A^) to DNS data with initial 
separation Aq = ijf/i (a) and Aq = rjf/2 (b). The DNS data is shown in dashed lines, black (forward) 
and red (backward). The model data is given in black for three values of fco: —0.85, —0.75, —0.65 (top to 
bottom). The lines are very close together, indicating a weak dependence of (A^) on fco- The diamonds 
(♦) demark the threshold time tc, which is the largest time for which one can expect agreement of the 
model with the data. 
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Figure 9: Numerically integrated mean square separation (A^) (solid black) compared to the simple 
approximation (56) (dashed black) and DNS data (solid red) for all 10 initial separations, using the best 
fit values r/t = 1.2 and fco = —0.75. 
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This is shown in fig. 9, along with the present model and the DNS data. It easily satisfies the two given 
limits (A^) = Ao for t — and Richardson's law (46) for t — » cx). In general, the present model and 
this simple approximation are more similar to each other than each compared with the DNS data. For 
the best fit cases with initial separation Aq G {rjf /A,rif /2} (3''* and 4**^ from the bottom in fig. 9), the 
present model is a better fit than the given approximation (56). 

Thus, this approximation correctly traces the qualitative features of the present model and hence, 
provides a useful analytical approximation which can be used as an estimate. However, its use for 
accurate data fitting remains questionable. 



3.4 Time asymmetry in two dimensional turbulence 




Figure 10: Ratio of mean square separations (Af^^^) / (A^^^) for all initial separations from the DNS 
data. The lines show the largest initial separation to the smallest through the color spectrum red, blue, 
green, yellow, orange. If color is not available, this corresponds roughly to bottom to top order at i « T^. 
The Batchelor time ts is indicated by triangles (a) and the threshold time tc by diamonds (♦). Note 
that the "best fit" cases with initial separation Aq G {'7//4, ?7//2} are two (green) curves near the highest 
achieved ratio. 

Looking at fig. 7, it seems that there is no significant difference in mean square separation between 
the forward and backward case in the investigated two dimensional stationary turbulence DNS. 

However, the logarithmic scale is deceiving, as the plot of the ratio (Aj^^) / (A^^^) in fig. 10 shows. 
It indicates that the the mean square separation in the forward case is up to 22% larger than the backward 
case. 

The time asymmetry for a 3D FTV"^ experiment with Re^ ~ 170 has been documented in the 
literature by giving the ratio of the Richardson constants 3/ = GA,fwd and gb = GA.bwd from a fit of 
(56) [3]. In this 3D case, the observed ratio was gb/g/ ~ 2.1, i.e. the backward case separated stronger 
than the forward case. This ratio has been confirmed in the same work by a DNS with Re^ ~ 280. 

To be able to draw a more direct comparison to the 2D case, we shall plot the same ratio in fig. 11, 
obtained from (56), over time: 

//A2 \l/3 a2/3\^ 



9b \ /A2 \^/^ _ A^/^ 



(57) 



This relation is useful, since it takes the qualitative dependence of Richardson's law on the initial sepa- 
ration Aq into account, as described in §3.3. Fig. 11 shows that this ratio is larger than unity for most 
initial separations and times. Furthermore, with the exception of the Ao = ?///16 dataset, around (and 
between) the critical times ts and tc, the ratio is between 1.05 and 1.25. Hence, the best estimate from 
fig. 10 and fig. 11 for the ratio of Richardson's constants in this 2D DNS is: 

^ = (1.15 ±0.10) . (58) 
9b 



^Particle Tracking Velocimetry 
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Figure 11: Ratio of Richardson's constants gf/gb from the approximation (56) for all initial separations 
from the DNS data. The lines show the largest initial separation to the smallest through the color 
spectrum red, blue, green, yellow, orange. If color is not available, this corresponds roughly to bottom 
to top order at t « T^. The Batchelor time is indicated by triangles (a) and the threshold time tc 
by diamonds (♦). 
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Note that in 2D the forward separation is stronger while the effect is opposite in the 3D experiment, where 
backward separation dominates. However, the effect in 2D with ~ 15% difference is not as prominant as 
in 3D 100%). The figures fig. 10 and fig. 11 suggest that the ratio is different from unity. However, 
given the uncertainty of the ratio, one cannot conclude this for certain from this kind of fit. 

3.5 Richardson's constant and improvement of time asymmetry fit 

While it was possible to estimate the ratio of Richardson's constants from the time scries of mean square 
separation ^A^(i)) for both the forward and backward case, it is well known that obtaining Richardson's 
constant itself from this kind of data is much more problematic [4, 18]. Alternative suggestions of ex- 
tracting Richardson's constant include statstics of doubling times [4] and investigation of mean diffusivity 
depending on mean separation [18]. 

We will follow the latter approach. Firstly, note that Richardon's constant Ga is usually defined in 
terms of the law: 

(A2(t)) =GAet\ (59) 
where e is the mean energy dissipation of the flow. From this, one can easily derive the relation 

^^{A\t))=3{GAef' {A'{t)f\ (60) 

which is valid for all times if one assumes that the Richardson law holds. This relation also agrees with 
(56), taking the qualitative dependence on initial separation into account. 

To analyze the DNS data for (A^(i)), the value of the time derivative is obtained by using a finite 
difference scheme of second order. This data is given in logarithmically equidistant bins and therefore, 
the scheme needs to account for variable distances between the data points: 



df 1 



dt dt_i + dti 



d<i ^ ^ V dh At^i) ^ ' dt_i ^ 



clt_i) 



0[di3] , (61) 



where /(t) is an arbitrary function, d<_i is the distance to the previous timestcp and dti is the distance 
to the next timestep. 

However, it was not possible to observe relation (60) in the DNS data as can be seen in fig. 12. 
Instead evidence for exponential growth is found. NicoUcau and Yu [18] found similar behavior for small 
separations in a 3D kinematic simulation with a large inertial range. While they also observed a region 
(X (A^)^^'^, we suspect that the absence of this observation in fig. 12 is due to the limited width of the 
inertial range. 
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Figure 12: DifFusivity d/dt(A^(t)) vs. mean square separation (A^(t)) for all initial separations from 
the forward DNS data (smallest to largest left to right). The expected slope of 2/3 as predicted by (60) 
cannot be observed. Instead the best fit (straight line) suggests the exponent 1 indicating exponential 
growth. 



Note that independently of the slope of fig. 12, the ratio of change of mean square separation also 
gives the ratio of Richardson's constants: 

where we used (59) in the second last equality. Taking finite initial separations into account by using 
(56) gives the same result for large enough times. 




Figure 13: Ratio of Richardson's constants gj/gb as given by (62) for all initial separations from the 
DNS data. The lines show the largest initial separation to the smallest through the color spectrum red, 
blue, green, yellow, orange. If color is not available, this corresponds roughly to bottom to top order at 
t « Tc- The Batchelor time Ib is indicated by triangles (A) and the threshold time tc by diamonds (♦). 

Plotting this ratio using the DNS data in fig. 13 looks remarkably like fig. 11 and is also consistent 
with the value of 

^ = (1.15 ±0.10) . (63) 
gb 

found earlier. Despite the fact that fig. 12 suggests non-Richardson exponential growth of (A^(t)) due to 
the limited power-law range of length scales, this still is evidence that the forward separation is happening 
faster than the backward separation. 
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3.5.1 Conditioning on separation scale 

NicoUeau and Yu [18] have demonstrated that it is possible to have exponential and Richardson like 
growth in the same flow. Our suspicion is that the Richardson like behavior is still present in the DNS 
flow, but overshadowed by contaminations from scales outside the inertial range, which are present in 
every average over all Lagrangian trajectories. 

Going back to Richardson's original concept that the diffusivity is purely scale dependent [20], we 
are going to subject (60) to conditioning on separation scales A: 

'a\ ^ 3 (Ga ef' {A'\Af' ^ 3 (Ga ef' A^'^ . (64) 



dt 

Note that this average is entirely different fr'om the previous ones. While in previous sections of this 
article, the averaging was carried out by sampling the set of Lagrangian trajectories for each timestep in 
the DNS (i.e. ensemble average, conditioned on DNS time t), it is now completely independent of time 
and instead averages over all particles pairs that have a certain separation A at any given time. Hence, 
temporal information is lost for averages of this kind, but the conditioning on scales ensures that scale 
contamination is completely eliminated. 

Since this approach docs not directly use time series data, the finite difference approach of obtaining 
the time derivative is not practical. Instead, let us consider the relative motion of particle 1 and 2 in a 
pair: 

A ^ rr^ ^ {xi - X2) ■ {ui - 112) -^ A- fRK\ 

— A = — ^{xi - X2Y = ^ = (Ml - U2) • A = MA , (65) 

where and Ui arc particle position and velocity and A = (xl — 3^2)/ A is the unit vector of particle 
separation, ma is then the longitudinal velocity increment'* of a particle pair, which can be calculated 
from the velocity field at every instance. Thus, we obtain the change of mean square separation from 

^)=2(Ama). (66) 

We can now without problem introduce conditioning on separation: 

'dA2 



dt 



A ) = 2 (Ama|A) = 2A (ijaIA) . (67) 



Note that the last equality holds for conditioning using infinitesimal A bins, but for the purposes of ob- 
taining data using finite width separation bins, more accuracy is obtained by using the second expression. 
Thus, combining (64) and (67) we arrive at a relation with can be tested against the DNS data: 

(AmaIA) = ^ (Gae)'/' A^/^ (68) 

The forward and backward DNS data is shown in fig. 14 and fig. 15. Both demonstrate the validity 
of (68) within the inertial range. It seems that this behavior is extended for separations even larger than 
C, while for A < 77/ the exponential separation typical of ballistic separation is observed. 

As expected, the relation (68) is independent of initial separation which manifests itself in a good 
collapse of all datasets, proving it a useful tool to circumvent the problem of dependence on finite initial 
separations. 

The best fits of Richardson's constants from (68) for both forward and backward separation are given 
in table 3. The overall estimate for Richardson's constant in this 2D turbulence is gf = (1.066 ± 0.020) 
for forward separation and gt ~ (0.999 ± 0.007) for backward separation. This leads to a ratio of 

^ = (1.07 ±0.03) (69) 
9b 



''Note that this longitudinal velocity increment in not the same as the one generally used in e.g. Kolmogoroff's 4/5 
law. The latter velocity increment is an average over all positions or ensembles at a particular time, whereas the velocity 
increment defined in (65) is an average over specific particle pairs at different times. A comparison of dimensions could lead 
to a suspicion that the two could have the same value. However, this is unlikely since Kolmogoroff's 4/5 law has opposite 
signs in 2D and 3D, while ma is positive in both 2D and 3D. 
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Figure 14: (Aua|A) vs. pair separation A for some initial separations from the forward DNS data. 
The expected slope of 4/3 as predicted by (68) can clearly be observed for at least the incrtial range 
(—2 « log 77/ < log A < log£ « —0.6) with gf = (1.066 ± 0.020). For smaller separations, the expected 
exponential growth for ballistic separation with slope oc A^ is also observed. Note that separations 
A < Aq have only a small number of pairs contributing to the statistics and therefore, exhibit a large 
scatter. 
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Figure 15: (Awa|A) vs. pair separation A for some initial separations from the backward DNS data. 
The expected slope of 4/3 as predicted by (68) can clearly be observed for at least the inertial range 
(—2 « log?7/ < log A < log£ « —0.6) with gt = (0.999 ± 0.007). For smaller separations, the expected 
exponential growth for ballistic separation with slope oc A^ is also observed. Note that separations 
A < Aq have only a small number of pairs contributing to the statistics and therefore, exhibit a large 
scatter. 



Ao 


9f 


9b 


9f/9b 


77/ /16 


1.061 ±0.004 


1.005 ±0.002 


1.06 ±0.01 




1.084 ±0.004 


0.996 ±0.003 


1.09 ±0.01 




1.080 ±0.004 


0.999 ±0.002 


1.08 ±0.01 


rjf/2 


1.053 ±0.004 


0.996 ±0.002 


1.06 ±0.01 




1.050 ±0.004 


1.001 ±0.003 


1.05 ±0.01 


combined 


1.066 ±0.020 


0.999 ±0.007 


1.07 ±0.03 



Table 3: Results from fit of (68) to the DNS data. Given are the best fit values of Richardson's constants 
for forward and backward DNS data and confidence intervals as found by gnuplot 4 . for separations 
within the inertial range (—2 « logTy/ < logA < logC « —0.6). Data with Aq > 77/ did not have 
good enough statistics within the inertial range to be considered for this fit. The uncertainty given for 
the individual fits is the asymptotic standard error given by gnuplot and the uncertainty given for the 
combined (averaged) values is the standard deviation (t„_i of the distribution of values above plus the 
average error above. 
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which is compatible with the earlier value (63). 

A further consistency check can be performed by fitting a constant to the ratio of Richardson constants 
as obtained from 



(A^aIA) 



fwd 



GA.bwd 



1/3 



1/3 



5/ ^ f{AuA\A)j^ 

m V(A«A|A)b^d 



(70) 




Figure 16: Ratio of Richardson's constants gf/gt from (70) for initial separations from the DNS data 
with Ao < rjf within the inertial range (—2 w log 77/ < log A < log£ w —0.6). The line colors are the 
same as in fig. 13. In the absence of color note that the ratio is > 1. The solid black line shows the best 
fit constant gj/gt = 1-10. 



Ao 


9f/9b 


77/ /16 


1.120 ±0.006 


Vf/8 


1.109 ±0.004 


W/4 


1.106 ±0.003 


W/2 


1.121 ±0.005 


If 


1.069 ±0.005 


combined 


1.105 ±0.022 



Table 4: Best fit values of ratios of Richardson's constants for forward and backward DNS data and 
confidence intervals as found by gmiplot 4.0 for separations within the inertial range (—2 w log?// < 
log A < log£ w —0.6). Data with Aq > ijf did not have good enough statistics within the inertial range 
to be considered for this fit. The uncertainty given for the individual fits is the asymptotic standard 
error given by gnuplot and the uncertainty given for the combined (averaged) values is the standard 
deviation cr„_i of the distribution of values above plus the average error above. 

The results from fitting (70) to the data are compatible with the previous value of gf/gt (see table 4). 
It can also be seen from fig. 16 that the ratio gf/gb is larger than unity for all five initial separations 
shown. We conclude that the ratio of Richardson's constants for this 2D DNS of isotropic homogeneous 
turbulence with resolution N = 3072 is the weighted average of the values given in table 3 and table 4: 

^ = (1.09 ±0.03). (71) 
9b 

4 Discussion 

The main limitation of the present model has been shown repeatedly in this article. It currently models 
an infinite inertial range with a multi-scale stagnation point topology of infinitely wide scale range. 
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However, in order for the model to accurately describe the systems available in experiment and DNS, it 
would be necessary to include finite range effects. 

Some suggestions are to incorporate the separation behavior similar to Brownian motion for scales 
larger than the integral scale (§3.1) and to include dissipative effects for separations that are smaller 
than the distance between stagnations points. 

All results obtained from comparison of DNS with the model need to be discussed in this light. Thus, 
the value obtained by GV04 of q;gv = 1.3 might turn out to be in agreement with Richardson's prediction 
of ctQY = 1.5, when one takes into account that the range of scales in the multi-scale stagnation point 
topology is finite (§2.1.5). 

These finite range effects are also the reason that Richardson's law is not clearly present in the DNS 
with a resolution of = 3072. However, this does not mean that one cannot interpret the DNS data at 
all. In §3.5.1, we were able to use a novel scale-dependent approach which confirms Richardson's scalings 
[20] in our finite range DNS even though the usual diagnostics and statistics do not show them clearly. 
Richardson's constants for the forward and backward case of this 2D DNS of isotropic homogeneous 
turbulence with resolution N = 3072 are found to be 

GA.fwd = (1.066 ± 0.020), (72) 
GA.bwd (0.999 ± 0.007) . (73) 

These values are approximately 4 times lower than the Richardson constant g « 3.8 for a similar 2D 
DNS with similar resolution as obtained by Boffetta and Sokolov [4] using exit time statistics and 
approximately twice as large as the value of g ~ 0.5 obtained from fitting the cx law to a 2D turbulence 
experiment by JuUien et al. [16]. The spread in obtained Richardson constants may be due to non- 
universality stemming from different characteristics of the investigated flows, but might also stem from 
the varying approaches in obtaining them. The notion of a non-universal Richardson constant is plausible 
and supported by 

Ga oc Cl Cl'"^ , (74) 

which can be derived from equations (6), (11) and (46). This predicts that the Richardson constant 
depends on the stagnation point number Cg and the constant Gb, both of which arc not necessarily 
universal. 

Finally, the investigation of the two dimensional DNS in §3.4 has shown qualitatively that the forward- 
/backward time asymmetry of turbulent pair separation is opposite in 2D to what it is in 3D. Furthermore, 
utilizing conditioning on scales rather than on time yields a sound quantitative value for the ratio of 
Richardson's constants in the present 2D turbulence DNS (§3.5.1): 

^ = (1.09 ±0.03) . (75) 
gb 

This can be compared to a value for 3D turbulence obtained experimentally in a previous study [3] : 

present 2D DNS : — = (0.92 ± 0.03) (76) 

9f 

3D experiment [3] : — = (2.1 ± 0.3) (77) 

91 

Berg et al. [3] explained the time asymmetry in terms of the positive sign of the mean second eigenvalue 
of the strain tensor in 3D turbulence. Their argument would imply that gf /gt = 1 in incompressible 2D 
turbulence, which we have shown not to be the case. However, this does not imply that their argument 
is wrong. It implies that their mechanism cannot be the only one contributing to the observed time 
asymmetry of pair separation. We shall attempt to give a qualitative explanation for what might be a 
different and perhaps additional mechanism. 

The mechanism causing this time asymmetry in two dimensional turbulence could have its origin 
in the asymmetry which exists between forward and inverse energy cascade. The latter is related to 
merging and thereby growing eddies in 2D turbulence, whereas the former is related to eddies breaking 
up into smaller eddies. 

It is known that in two dimensional turbulence, small vortices within each others proximity can merge 
into a vortex of a larger scale [7]. This idea of merging vortices has been explored for many decades 
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Figure 17: Diagram illustrating the merging of eddies. Shown on the left are two neighboring eddies 
with their elliptic (■) and their hyperbolic (•) stagnation points and a characteristic streamline. The 
right shows the larger eddy resulting from this merging process and its elliptic stagnation point. The 
merging involves the transformation of two elliptic and one hyperbolic to one elliptic stagnation point. 

[1, 12]. Picking up the notion of a particle pair's patron from §2, let us point out that it might be 
possible for a particle pair to move to a larger patron without encountering a hyperbolic stagnation 
point as presented in an earlier section. This alternative merging process would involve the annihilation 
of the small vortices and therefore the disappearance of their hyperbolic and elliptic zero-acceleration 
points and the creation of an elliptic zero-acceleration point for the emerging larger scale vortex, shown 
in fig. 17. This dynamic process of a particle pair moving to a larger patron cannot be explicitly included 
in our model which is based on particle pairs encountering hyperbolic stagnation points. However, it is 
possible to include the effect of this process in an effective value of either of the constants Cs, Cb or both, 
which will be smaller if pairs separate whilst eddies break up into smaller eddies (as in our time-reversed 
2D turbulence). Conversely, if pairs separate whilst eddies merge into larger eddies (as in forward time 
inverse-cascading 2D turbulence - see fig. 17), the value of either Cs and/or Cb would be larger. 

Summarizing, apart from the burst-like processes of separation and convergence described in §2, there 
might be an additional mechanism of merging vortices that also contributes to a particle pair's separation 
in 2D turbulence. This additional route to larger separation would strengthen the separation process in 
the natural time direction (forward), and would be reversed and therefore weakening the separation in 
the backwards case. Note that this process would be happening on time scales larger than the life time 
of the involved stagnation points. Note also that this process has not (yet) been directly observed in 
the context of pair separation and thus, is merely a suggested notion which consistently ties together all 
observations. 

The inverse energy cascade of two dimensional turbulence is understood to describe the same phe- 
nomenon [7] of smaller vortices merging into larger ones. While nothing more than a handwaving 
argument, transferring this notion to three dimensional turbulence, where the energy cascade is associ- 
ated with larger vortices breaking up into smaller ones, the separation process is depleted in the natural 
progression of time and strengthend in the backward case, giving a consistent qualitative picture for the 
observed time asymmetry in 2D and 3D turbulence. 

5 Conclusions 

An improved model of particle pair separation in flows with a multi scale stagnation point topology was 
devised, thus setting its predecessor [14] on a sound mathematical foundation. The new model is able to 
model the observed time asymmetry provided that the effects of vortex- merging (fig. 17) can be taken 
into account in effective values of and/or Cb- Furthermore, it has been argued that correction terms 
of higher order derivatives are necessary in the PDF evolution equation (14) when considering separation 
moments of third order or higher. 

The limitations and possible improvements to the present model are discussed, and the Richardson 
constant for the isotropic homogeneous 2D DNS turbulence is found to be of the same order of magnitude 
as previous comparable values [4, 16]. 

The time asymmetry is found to have opposing effects in two and three dimensional turbulence. It is 
suggested that the stronger forward separation in 2D turbulence might be caused by merging of eddies. 
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Furthermore, the observed asymmetry is consistent with the assumption that the direction of the energy 
cascade in 2D and 3D turbulence is directly correlated with the direction of the time asymmetry. 
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